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Abstract 

A mechanical model of swimming and flying in an incompressible viscous fluid in the absence 
of gravity is studied on the basis of assumed equations of motion. The system is modeled as an 
assembly of rigid spheres subject to elastic direct interactions and to periodic actuating forces 
which sum to zero. Hydrodynamic interactions are taken into account in the virtual mass matrix 
and in the friction matrix of the assembly. An equation of motion is derived for the velocity of the 
geometric center of the assembly. The mean power is calculated as the mean rate of dissipation. 
The full range of viscosity is covered, so that the theory can be applied to the flying of birds, as 
well as to the swimming of fish or bacteria. As an example a system of three equal spheres moving 
along a common axis is studied. 
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I. INTRODUCTION 


The swimming of fish and the flying of birds continue to pose challenging theoretical 
problems. The physics of bird flight was hrst studied in detail by Otto Lilienthal in the 
nineteenth century 1^. Since then, signihcant progress has been made in many years of 
dedicated research [21]-[51]. 

The goal of theory is to calculate the time-averaged speed and power for given periodic 
shape variations of the body, at least for a simple model system. It is assumed that the 
motion of the fluid is well described by the Navier-Stokes equations for an incompressible 
viscous fluid. On average over a period the force exerted by the body on the fluid vanishes, 
so that thrust and drag cancel. In early work by Lighthill [^ and Wu [^ the thrust and 
power were calculated approximately as functions of the speed on the basis of potential flow 
theory for a planar st™. This work and subsequent developments have been reviewed by 
Childress [^, by Wu |8f,[^, and by Sparenberg [l^. However, an independent calculation 
of the mean speed for given periodic shape variations is still lacking. Measurement of the 
power consumption has led to a surprisingly small friction coefficient, much smaller than 
that of an inert body, as was hrst observed by Gray flU ]. 

It was hrst shown by Taylor [l^ that in the slow swimming of a microorganism the 
calculation of thrust can be circumvented. In this limiting case one can use the time- 
independent Stokes equations. The mean swimming velocity and mean rate of dissipation 


then follow from a purely kinematic calculation [I3|,[l^. For small amplitude swimming 
both quantities are quadratic in the amplitude of the stroke to lowest order. For a simple 
system, where the body is modeled as an assembly of rigid spheres held together by direct 
interaction forces and subject to periodic actuating forces which sum to zero, we have shown 
that in the high viscosity limit the swimming velocity and power can be calculated for any 
amplitude of stroke from kinematics alone 15|].(l^. 

In the following we investigate questions of thrust, velocity, and power for swimming or 
flying in a fluid of any viscosity, including the limit of low viscosity, for the same mechanical 
model as before. We assume for simplicity that the spheres experience Stokes friction. 
In addition we incorporate hydrodynamic interactions via virtual mass effects, as found 
from potential flow theory. We use Hamilton’s equations of motion with added damping 
terms. In the limit of high viscosity, where resistive forces dominate, the earlier results 
are recovered. The model provides valuable insight also in the limit of low viscosity, where 
reactive forces dominate. In that regime the motion is dominated by virtual mass effects. 
Bernoulli forces and modihed linear friction should be taken into account in a more realistic 
model. Nonetheless, the principle of the calculation, which exploits elimination of the fluid 
degrees of freedom, remains valid. 

The flow is assumed to be laminar at all times. It is now realized that the boundary layer 
of swimming fish is laminar even at high Reynolds number [^. Virtual mass effects were 
discussed earlier by Lighthill 11111 . The numerical modeling of animal swimming and flight 
was reviewed by Deng et ah |l8j . 

As an example a system of three equal spheres moving along a common axis is studied. 
For this simple system the mean swimming speed and mean power to second order in the 
amplitude of stroke can be evaluated analytically. The solution to a corresponding eigenvalue 
problern provides the optimal stroke to this order, as we found elsewhere in the resistive 
regime [l5|. 

In our model the mean thrust, i.e. the frictional force exerted on the fluid averaged 
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over a period in periodic swimming, vanishes identically. We find that the velocity of the 
geometric center of the assembly is driven by a different force, which we call the impetus. It 
has both a reactive and a resistive component. The impetus determines the center velocity 
with retardation. The mean impetus does not vanish. 

It is known for small amplitude swimming in the resistive regime that the mean power 
is directly proportional to the mean velocity. We find for our example that the relation 
between mean power and mean velocity is nearly linear also for large amplitude swimming. 
Presumably the near linearity holds also for other systems in the whole regime of viscosity. 
If true, this would resolve the so-called Gray paradox [^, which is based on the mistaken 
notion that the power is quadratic in the velocity, as in Stokes friction. 


II. EQUATIONS OF MOTION 


We consider a set of N rigid spheres of radii oi,..., oat and masses m^i, ...^rUpN, centered 
at positions R = {Ri,R n), and immersed in an incompressible viscous fluid of shear 
viscosity rj and mass density p. The fluid is of infinite extent in all directions. The flow 
velocity v{r,t) and pressure p{r,t) of the fluid are assumed to satisfy the Navier-Stokes 
equations 

p [— -h V ■ Vr;] = — Vp, V ■ t; = 0. (2.1) 

The flow velocity v is assumed to satisfy the no-slip boundary condition on the surface of 
the spheres. The fluid is set in motion by time-dependent motions of the spheres. At each 
time t the velocity field v{r, t) tends to zero at infinity, and the pressure p{r, t) tends to the 
constant ambient pressure po- 

As the spheres move in the fluid they experience a frictional force. In addition there may 
be applied forces E(f) = {Ei(t),EN(t)) and direct interaction forces which depend on 
the relative positions {Rj — Rk} of sphere centers. We shall assume that the sum of applied 
forces vanishes, so that 

N 

( 2 . 2 ) 


The sum of direct interaction forces vanishes owing to Newton’s third law. We assume that 
the frictional forces are linear in the sphere velocities, as given by low Reynolds number 
hydrodynamics on the slow timescale 19f|. The spheres are freely rotating so that there are 
no frictional torques. 

The forces exerted by pressure gradients resist instantaneous acceleration and give rise 


to virtual mass effects [201]. For a single sphere j immersed in infinite fluid the virtual mass 


would be where rufj = Anpa^/S is the mass of fluid displaced by the sphere. The 

virtual mass effect for a collection of N spheres is embodied in a {3N x 3A)-dimensional mass 
matrix m. This can be derived from potential flow theory by considering the irrotational 
flow pattern generated instantaneously by a set of sudden impulses S = (Si, ..., Sn) from a 
state of rest. The total corresponding kinetic energy is 


JC 


U • m ■ U, 


(2.3) 


where U = (C/i,..., U n) is the set of sphere velocities. The kinetic energy is a sum of two 
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contributions, /Cp + /C/, with 


^ E “wt'i ■ >^f = ^P f *. (2-4) 

2,.i 2 7v, 

where the integration is over the part of space occupied by fluid and 0(r) is the velocity 
potential corresponding to the set of instantaneous velocities U. The potential 0(r) is linear 
in the sphere velocities U = (C/i,C/ tv), so that the kinetic energy /Cj is a quadratic form 
in U. The contribution /Cj to the total kinetic energy dehnes the virtual mass. This depends 
parametrically on the positions R, leading to corresponding hydrodynamic interactions. 

Besides the mass matrix m it will be useful to dehne also the mobility matrix /x, the 
friction matrix and the inverse mass matrix w, 

H w = (2.5) 

The four matrices are symmetric and depend on the relative positions of the sphere centers. 
The sphere momenta, including the virtual mass contribution, are given by 

p = m ■ U, U = w ■ p. (2.6) 


Correspondingly we postulate the equations of motion 


dR dp die dVint ^ 

— = U, — =-C-U-+ E 

dt dt dR (9R 


(2.7) 


where /C is given by /C = |p ■ w ■ p and Vint is the potential of direct interaction forces. The 
partial derivative d/dR is taken at constant momenta p. It is clear that the equations of 
motion (2.7) have only limited validity. The frictional forces are assumed to be linear in the 
sphere velocities, and Basset memory forces are neglected. Nonetheless it is of interest to 
explore the consequences of the equations as they stand. 

We note that it follows from Eq. (2.7) that total dressed sphere momentum P = SjPj 
is not conserved but changes due to friction with the fluid. If an impulse S is imparted at 
time t = 0 to the spheres in a state of rest due to applied forces E = SS(t), then part of 
the momentum is transferred instantaneously to the fluid, reducing the sphere velocities to 
U(0+) = w ■ S, as can be seen by integration over an inhnitesimal time interval about t = 0 
and use of Eq. (2.6). The remaining momentum of the spheres is transferred to the fluid in 
the course of time by friction. The total momentum of spheres and fluid is conserved at all 
times. 


III. IMPETUS, CENTER VELOCITY, AND RATE OF DISSIPATION 

We are looking for a solution of the equations of motion where the center of the assembly 
moves on average with uniform translational velocity and the individual spheres perform 
periodic motions about the moving center. The velocity U{t) of the geometric center is 
dehned by 

1 

= — U ■ u„, a = {x, y, z) (3.1) 
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where the symbol Ux denotes a 3iV-dimensional vector with 1 on the x positions, 0 on the 
y,z positions, and cyclic. The 3iV-dimensional displacement vector d(t) = SN(t)) 

is dehned by 

Rj{t) = Rjo + f U{t') dt' + 5j{t)^ j = 1, (3-2) 

Jo 


with the property 


d(t) ■ u„ = 0, a = {x,y,z), 


(3.3) 


and with a set of equilibrium positions for which the direct interaction forces vanish. 

Correspondingly the velocity vector is decomposed as 


U = UisUy + d. 


(3.4) 


Summation over repeated greek indices is implied. By substitution into Eqs. (2.6) and (2.7) 
we hnd 

— [m ■ {UyUy + d)] + — + C ■ {UyUy + d) + = E(t). (3.5) 

Multiplying from the left with the vector Uq, we obtain 

— + — (uo ■ m ■ d) + + Uq, ■ ^ ■ d = 0, (3.6) 

with mass tensor AT and friction tensor Z dehned by 


= Ua ■ m ■ U/ 3 , Zafi = u„ ■ C ■ u^. (3.7) 

We have used the fact that m depends only on relative coordinates, so that Uq, ■ dfC/dR = 0. 
Also we have used Newton’s third law and Eq. (2.2). We note that in Eq. (3.6) the center 
velocity Uy occurs only in the two places explicitly shown. 

We rewrite Eq. (3.6) as 

^(A4-C/)+Z-C/ = X, (3.8) 

(Jjb 

with time-dependent impetus I{t) given by 

Zait) = (ua ■ m ■ d) - Uo ■ C ■ d. (3.9) 

The mean impetus, averaged over a period r, 

la = — [ Za{t) dt = —Uq, ■ C ■ d, (3.10) 

Jo 


does not vanish, even though the total mean force exerted by the huid on the assembly does 
vanish. The drag exerted on the assembly by the huid is 


Rex LIq ■ ^ ■ U ZafjUp Uq ■ ^ ■ d. 


(3.11) 


The thrust T is equal and opposite to the drag, T = —D. There are frictional forces 
on the spheres, but on time average the total drag vanishes, D = 0, T = 0, as follows 
from Eqs. (2.2), (2.7) and Newton’s third law for the interaction forces. We use the fact 
that the integral of Uq ■ dp/dt over a period vanishes by periodicity, as well as the property 
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Wa ■ dfC/dR = 0, mentioned below Eq. 
can also be expressed as 


(3.7). From D = Q it follows that the mean impetus 


2^0: ^OL^Up. 


(3.12) 

We may regard Eq. (3.8) as a balance between central and internal motion. The equation 
must be solved for the center velocity U{t) for given impetus X(t), which can be calculated 
from the displacement vector d{t). In the resistive limit the total drag vanishes at any time 
and the reactive forces vanish, so that then Ua = MapXp = —^ ■ d. 


21 


The displacement vector d(f) may be calculated from displacements in relative space 
by using a transformation to center and relative coordinates. The geometric center of the 
assembly is given by 

1 ^ 1 

Rj = — e„Ua ■ R (3.13) 


R = 


N 


i=i 


N 


with Cartesian unit vectors e„. We define relative coordinates {vj} with j = 1,..., iV — 1 as 

Vi = R 2 — Rl, r2 = i?3 — i?2; •••; 

vn-i = Rn — Rn-1, (3-14) 


and denote the corresponding {3N — 3)-vector r = (ri,..., rAr_i). The 3Wvector {R,r) is 
related to the vector R by a transformation matrix T according to 


(i?,r) = T-R (3.15) 

with explicit form given by Eqs. (3.13) and (3.14). The displacement vector d is derived 
from the displacement ^ in relative space as 

d = T-i-(0,^). (3.16) 


Therefore the impetus is determined by displacements in relative space. 

The time-dependent rate of dissipation can be expressed in the same matrix formalism. 
The rate of dissipation is given by 

P = U-C-U. (3.17) 

Once the center velocity U{t) has been calculated for known displacements we can also 
calculate the time-dependent rate of dissipation 'D{t) by use of Eq. (3.4). However, we 
derive an alternative expression which will be useful in the following. We can solve for the 
product C ■ U from the equation of motion 


|(-U)+C-U = F. 

where F is the vector of mechanical forces acting on the spheres, 

die dVi^, 
dR dR ^ ’ 


(3.18) 

(3.19) 


which has the property Uq ■ F = 0. Using this property we find for the rate of dissipation 


P=U.F-U.4(„,.U)=d.F-U.4(n,.U). 

Using here Eq. (3.18) again we can rewrite this as 


(3.20) 


with friction vector f„ 
limit 
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V = d ■ C ■ d + Uad ■ Uaiia • ^ (ui • U), (3.21) 

C ■ Uq. This generalizes an expression derived earlier in the resistive 
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IV. SMALL AMPLITUDE SWIMMING 


For vanishing displacements Eq. (3.8) has the solntion U = 0. By formal series expansion 
in powers of the displacement vector d we obtain a corresponding expansion of the center 
velocity 

= + (4.1) 

The hrst order velocity satishes the eqnation 

+ Z°i,uf> = -u„ . m“ . a - u„ . C“ ■ d, (4.2) 

where the snperscript 0 indicates that the qnantity is calcnlated for the conhgnration Rq. In 
particnlar for oscillating displacement, in complex notation 

d(t) = Re [d^e-'^*], (4.3) 

we hnd correspondingly 

Ui]} = [ - iujM^ + • (a;2m° + iuC°) ■ d^. (4.4) 

In this sitnation the hrst order velocity oscillates in time and vanishes on time average. 

Multiplying Eq. (3.8) by the mobility tensor M = Z~^ and expanding to second order 
we obtain a more complicated eqnation for the second order velocity It snfhces to 

derive an expression for the time-average over a period r = ‘Injuj^ 

C/(2) = i/’ u^‘^\t)dt. (4.5) 

Jo 

Using periodicity and the hrst order eqnation Eq. (4.2) we obtain the expression 

uF = . C'‘> ■ d + (4.6) 

Alternatively the expression can be derived directly from Eq. (3.12) by nse of = 
-Z°M(^)Z°. Here we can nse 

C^^^ = d-VClo, M(^) = d-VM|Q, (4.7) 

where V indicates the gradient operator in conhgnration space, and the notation indicates 
that the valne of the gradient is taken at Rq. The time average in the hrst term in Eq. (4.6) 
can be expressed as 

U/3 ■ ^ ■ dj, (4.8) 

with derivative friction matrix 

= Vf^, f^ = C ■ u^, (4.9) 

as introdnced earlier (l^ . 

In the second term in Eq. (4.6) we nse the identity 

ZayM^I3 = 5al3 (4T0) 


7 








to show that 


= —McfygjMsiS 


(4.11) 


with gradient vectors 

g? = ■ u^. (4.12) 

The hrst order velocity is eliminated by use of Eq. (4.4). Then the second order mean 
swimming velocity can be expressed as 

tiF=lRe[swd^V“(w)|j,.dJ, (4.13) 

with frequency-dependent matrix V“(a;) given by 

V“(a;) = Mo^^D^u), (4.14) 

with reduced derivative friction matrix 


with admittance tensor 
and impedance vector 


Df^{uj) = D^-g%s{ooMoo), 

Y{u:) = [-iujM + Z]~\ 


f<5(i^) = {-iojm + C) • Us. 
The matrix D^(ci;) has the property 

□^(cj) ■ Ua = 0. 


(4.15) 

(4.16) 

(4.17) 

(4.18) 


Since C depends only on relative coordinates = 0, and hence 

u„ • = 0, u„-g^ = 0. (4.19) 

As a consequence 

u„-V^(a;) = 0, V“(a;) ■ u ;3 = 0. (4.20) 

These properties generalize those derived earlier at zero frequency [l^ . 

To second order in the displacements the last term in Eq. (3.21) can be written as 

UaUa ■ ^(m • U) ^ ■ m° ■ d. (4.21) 

In the average over a period the hrst term on the right does not contribute. Hence for 
periodic motion the mean second order rate of dissipation can be expressed as 


•D(2) = 


Substituting Eq. (4.4) we hnd 


d • C° ■ d + Ui^'^d ■ f° - Ui^^Uo, ■ m° ■ d 


1 


dt. 


P(2) = - a;2 Re [d: • P(n;) ■ d^], 

with the complex matrix 

PM = C°->;°;.(^)f°(^)f». 

The matrix is symmetric and has the properties 

Ua ■ Re[P(a;)] = 0, Re[P(a;)] ■ = 0. 


(4.22) 


(4.23) 

(4.24) 

(4.25) 


The properties Eq. (4.20) and (4.25) allow us to reduce the dimension of the matrix de¬ 
scription by three by introducing center and relative coordinates. 
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V. VELOCITY MATRICES AND POWER MATRIX 


The transformation given by Eqs. (3.13) and (3.14) can be used to reduce the calculation 
of the second order swimming velocity and rate of dissipation to one in relative space. The 
matrices V“(a;) and P(ci;) are transformed to 

V“(n;) = T-V"(n;) •T-\ Pr(n;) = T ■ P(n;) • T’h (5.1) 

The hrst three rows of T consist of u^/iV and the hrst three columns of T“^ consist of Uq. 
It follows from the properties Eq. (4.20) and (4.25) that the first three rows and columns of 
the transformed matrices V^(a;) and Pt(i^) vanish identically. Hence in this representation 
we can drop the center coordinates and truncate the matrices by erasing the hrst three rows 
and columns. We denote the truncated (3iV — 3) x (3iV — 3)-matrices as \/^{u) and Pr(i^) 
and dehne displacements in relative space by 

(0,^J=T-d^. (5.2) 

With this notation the mean second order swimming velocity and rate of dissipation are 
given by 

IF = 1 Re ■ Ct • V?:(w). C„, 

^ = i Re C • Ct ■ Pt(w) • C. (5.3) 

with the matrix 

Ct= [T-i-T-i]: (5.4) 

This (3iV —3) X (3iV —3) dimensional matrix consists of numerical coefficients and is obtained 
from the corresponding 3N x 3N matrix by truncation, as indicated by the hnal hat symbol. 

We rewrite the expressions in Eq. (5.3) in a more convenient form with vectors and 
matrices consisting of dimensionless numbers. We introduce the complex dimensionless 
vector 

= I C. (5.5) 

where 6 is a typical length scale, and dehne 

A ^ ^ [(Ct ■ Pr(w)|()) +'i(CT-Pr(t^)lo) ]. (5.6) 

where the superscript s indicates the symmetric part, the superscript a the antisymmetric 
part, the single prime the real part, and the double prime the imaginary part. With the 
scalar product 

N-l 

(€V) = E ■ F (5.7) 

i=i 

the mean swimming velocity and mean rate of dissipation can then be expressed as 

kF= ia,KriB“ir). ^=l^v(riAir). (s.s) 
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The matrices B“ and A are hermitian. We call B“ the velocity matrix and A the power 
matrix. 

We ask for the stroke with maximum swimming velocity in a class of strokes with equal 
rate of dissipation for hxed values of the geometric parameters, hxed frequency w, and given 

values of viscosity r] and mass density p. Maximizing the quadratic form oob'^Ua'^ — 
with Lagrange multiplier A" we obtain the generalized eigenvalue problem 

= A“A^L (5.9) 

The eigenvalues {A"} are real. The maximum efficiency for motion in 

direction a is given by the maximum eigenvalue as 

(5.10) 

The set depends on the choice of Cartesian coordinate system. Fur¬ 

ther optimization may be possible by a rotation of axes. In particular cases a natural choice 
of axes will suggest itself. 


VI. POWER AND DISSIPATION 

We view the swimmer or flyer as a dynamical system in periodic motion, driven by 
actuating forces E(t) satisfying ■ E = 0 and E(t -|- r) = E(t). In an expansion in powers of 
the actuating forces the hrst order displacements are given by Eq. (4.3). To second order the 
mean swimming velocity and mean rate of dissipation are given by Eq. (5.8). In this section 
we relate the mean power supplied by the actuating forces to the mean rate of dissipation. 

The instantaneous power supplied by the actuating forces is given by 


P{t) = E{t) ■ U(t). 

(6.1) 

From Eq. (2.7) 


^(/C + I4„t) = -u-c-u + E-u. 

(6.2) 

Hence we hnd for periodic motion 


U-C-U = E-U. 

(6.3) 


This shows that on time average over a period the power is fully dissipated by friction. 

Since the mean thrust vanishes the usual dehnition of energy wastage makes no sense 
here. Instead we dehne the energy wastage S{t) as the difference 

£ = P-U -I. (6.4) 

In the theory of hsh swimming and bird flight the energy wastage has been associated with 
energy being lost to the vortical wake [^. We dehne the corresponding Eroude efficiency as 

VF = (6.5) 

This concept may be useful for the comparison of different strokes. 
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VII. HYDRODYNAMIC INTERACTIONS 


In order to perform explicit calculations we must specify the form of the hydrodynamic 
interactions appearing in the friction matrix and the mass matrix. In practice one uses 
approximate expressions which are presumed to be reasonably accurate in the range of 
distances considered. 

The friction matrix can be calculated from an approximation to the mobility matrix based 
on Oseen’s pair interaction 


19| . In this approximation the pair mobility tensor for the pair 


(j, fc) is given by 




67ir]aj 


16jk + 


Sttt] 


\^j 


Rk 


+ 


{Rj — Rk){Rj — Rk)) 
\Ri - Rk\^ 


(1 — Sjk)- 


(7.1) 


The mobility matrix ^ is composed of pair tensors, and the friction matrix ^ is its inverse. 

The calculation of the mass matrix m is based on potential flow theory. A dipole approx¬ 
imation to the mass matrix can be evaluated on the basis of an expression for the force on a 
sphere in a uniform flow in potential flow theory, as given by Landau and Lifshitz and 
by Batchelor [^ . 

In potential flow theory the flow velocity is expressed as r; = — V0 with a scalar potential 
(f), which satishes Laplace’s equation = 0 by incompressibility. A sphere of radius a, 
centered at the origin and moving with velocity [/ in a fluid at rest generates a potential 

^ -U, r> a, (7.2) 


corresponding to the dipole moment 


9r/ = 2 ® 


If the sphere is placed in a uniform flow Vo this is modihed to [24 


(7.3) 


9 = 5 - Do). (7.4) 

For a collection of N spheres in a fluid at rest at inhnity the velocities and dipole moments 
are related in dipole approximation as 

Uj = ^ qj + ^^Fjk-Qk, j = (7.5) 

with dipole interaction tensor 

_1 _|_ 3?^^ 

F,k = F{R,-Rk), Fir) = --, (7.6) 

where r = r/r. We abbreviate Eq. (7.5) as 

U = A-^-q, q = A-U. (7.7) 

The velocity of sphere j after a sudden impulse Sj from a state of rest is given by [^.[^ 
1 3 

{rripj + -TUfj) Uj = Sj + -rrifj ^ Fjk ■ Qk, J = 1, N. (7.8) 
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Substituting from Eq. (7.7) and solving for the velocities we find by use of Eq. (7.5) and 
mfj = 4:71 pa^/3 


U = w ■ s 

(7.9) 

with matrix 


w = [rUp — m/ + 47 rp.A] , 

(7.10) 

where the matrices trip and are diagonal with elements mpjl and mjjl. 

effective mass matrix is 

The approximate 

m = trip — ru/ + AirpA. 

(7.11) 

If the spheres are neutrally buoyant one has simply 


m = 471 p A., w = 

47ip 

(7.12) 


In our application we shall consider this case. 


VIII. THREE-SPHERE SWIMMER OR FLYER 


The simplest application of the theory is to a three-sphere swimmer or flyer with the three 
spheres aligned on the x axis, as studied by Golestanian and Ajdari jUf in the resistive limit. 
The spheres move along the x axis, and the y and z coordinates can be ignored jl^. There 
are only two relative coordinates ri = X 2 — X 1 and r 2 = X 3 — X 2 , and the relevant parts of the 
matrices and A are two-dimensional. The relevant parts are denoted as and A^. In 
the bilinear theory we consider a point ro in r-space with coordinates (di, ^ 2 ), corresponding 
to the conhguration Rq of the rest system. As an example we consider the case of equal-sized 
spheres with oi = 02 = 03 = a, equal masses mi = m 2 = m 3 = m = dvrpa^/S, and equal 
distances between centers di = ^2 = d. 

For this case analytic expressions for the matrices B** and A^ can be derived, but are 
too complicated to be presented. In the high viscosity limit the expressions reduce to those 
derived previously [l5|. The matrices, dehned with b = a in Eq. (5.6), depend only on the 
ratio d/a and the dimensionless viscosity 
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V* = — 


( 8 . 1 ) 


The two eigenvalues A± = ±A+, as well as the corresponding eigenvectors = (1,^±) with 
also depend only on these two variables. The dependence on the viscosity is 
surprisingly slight over the whole range of values. The absolute value |^+| is close to unity 
over the whole range. In Fig. 1 we show the variation of A+ and |^+| with p^ for d = 5a. 
The argument of increases slightly from 0.6278 tt at = 0 to 0.6285 tt at = 10®. 

In the bilinear theory the optimal orbit (ri(t),r 2 (t)) in relative space is given by r{t) = 
Vq + ^o(^) with Vq = (d, d) and 


^o(^) = Re exp(—icut)]. 


( 8 . 2 ) 


with amplitude factor e. 
given by 


do(f) 


The corresponding displacement vector in 





conhguration space is 


(8.3) 
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In Fig. 2 we show a snapshot of the spheres and their velocities in the instantaneous rest 
frame at f = 0 for £ = 3, d = 5a and r]^ = 0.01. In Fig. 1 of Ref. 15 we showed the elliptical 
orbit in relative space for d = 5a and e = 0.1 in the limit of high viscosity. In Fig. 3 we 
compare this with the corresponding orbit in the limit of zero viscosity, corresponding to 
small r;*. The two plots are indistinguishable on the scale of the hgure. 

From the periodic displacement do(t) we can calculate the instantaneous swimming ve¬ 
locity U{t) as a series of harmonics from Eq. (3.8). The zeroth harmonic yields the mean 
swimming velocity U. From U{t) and do(t) we can calculate the time-dependent rate of 
dissipation T>{t), and hence the mean V, by use of Eq. (3.17). 

In Fig. 4 we show the reduced mean swimming velocity U/{e‘^ua) as a function of e for 
d = 5a and r/* = 0.01. In Fig. 5 we show the reduced mean power P/{e‘^r]u‘^a^) vs. the 
reduced mean swimming velocity U/{e^uja) in the range 0 < e < 3. In Fig. 6 we show the 
efficiency Eqp = rjojo^u/P as a function of £. The efficiency increases monotonically with 
the amplitude factor. In Fig. 7 we show the time-dependence of the impetus X{t) and the 
center velocity U{t) as functions of time during a period for e: = 3. We also plot separately 
the resistive contribution to the impetus. At this value of p* this is much smaller than the 
reactive part. It is noteworthy that the variations in time of the center velocity are much 
smaller than those of the impetus. The center velocity follows the impetus with an after¬ 
effect. In Fig. 8 we show the absolute value of the Fourier coefficients /„ of the velocity 
7/(t), normalized to /o = 1, for e = 3. This shows that only a small number of harmonics 
contribute appreciably. 

The equality of mean power and mean rate of dissipation given by Eq. (6.3) can be 
checked numerically. The Froude efficiency rjF, dehned in Eq. (6.5), for these three values 
of e is 0.0004, 0.0020, 0.0057, respectively. 

It is of interest to compare the above results with values obtained from the numerical 
solution of the equations of motion Eq. (2.7) with hydrodynamic interactions given by Eqs. 
(7.1) and (7.12) and with prescribed oscillating actuating forces. In order to stabilize the 
system we consider microswimmers with internal harmonic interactions 2^. In matrix form 
the forces may be expressed as 


F = H • (R - Ro) + E, 


(8.4) 


with a real symmetric matrix H with the property H ■ = 0. The actuating forces {Ej{t)} 

can be chosen to correspond to the eigenvector with maximum eigenvalue in the problem 
Eq. (5.9). The hrst term in Eq. (8.4) represents a harmonic approximation to the direct 
interactions. We use harmonic interactions given by the 3 x 3-matrix 


(8.5) 


with elastic constant k. This corresponds to nearest neighbor interactions of equal strength k 
between the three spheres. The stiffness of the assembly is characterized by the dimensionless 
number a dehned by 

k 

a= -. ( 8 . 6 ) 



TTriau 


We use the hrst order equations of motion. 


dt 


= U«, 


m 


dUW 

dt 


= -C° ■ + H ■ R(^) + Eo, 


(8.7) 
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to calculate the actuating forces Eo(t) corresponding to the optimal linear motion given by 
Eqs. (8.3) and (4.4). These have the property Uq, ■ Eq = 0, so that the sum of actuating 
forces vanishes. We choose initial conditions corresponding to the rest conhguration 

a;i( 0 ) = 0 , 3:2(0) = d, 3:3(0) = 2d, 

f^i(O) = 0, 172(0) = 0, 173(0) = 0. (8.8) 

The kinetic energy term in Eq. (2.7) makes the direct numerical solution of the equations 
of motion time-consuming. Instead we use an iterative procedure, neglecting the kinetic 
energy term in Erst approximation. Then the equations can be solved by a fast procedure. 
From the solution the kinetic forces, defined as —dJC/dR, can be calculated as a function of 
time. In the next step we include the kinetic forces and repeat the procedure. The kinetic 
forces are small compared to the actuating forces and the solution converges after a few 
iterations. 

In Fig. 9 we show the numerical solution of the equations of motion Eqs. (2.6) and 
(2.7) with forces given by Eq. (8.4) for d = 5a, viscosity 77 * = 0.01, stiffness a = 10, and 
amplitude factor e = 1.5 for the first fifty periods. In Fig. 10 we show the mean value of the 
kinetic energy for successive periods. The orbit for the last period hardly differs from the 
ellipse given by Eq. (8.3), as shown in Fig. 11. The mean swimming velocity and the mean 
rate of dissipation can be calculated as time-averages over the last period. The efficiency is 
Et = 94 X 10“^, equal to the value calculated from the periodic orbit by use of Eq. (3.8) 
for displacements do(t) with e = 1.5, shown in Fig. 6 . 


IX. DISCUSSION 

In general the performance of a swimmer or flyer can be measured in terms of the dimen¬ 
sionless efficiency Et, defined as the ratio 

rjua^U 


where a is a conveniently chosen length scale, U is the mean speed and P is the mean 
power, averaged over a period r = 2ti/uj. We note that the lower bound of the efficiency 
Et vanishes, since for a periodic motion in relative phase space which is time-reversible the 
mean swimming velocity vanishes. A striking result of the present analysis is that for the 
simple three-sphere model, for which analytic calculations can be performed, the maximum 
efficiency ETmax is nearly independent of the dimensionless viscosity 77 * = rj/{uja?p), as shown 
in Fig. 1, see Eq. (5.10). As a consequence the mean speed is nearly inversely proportional 
to the shear viscosity 77 for given power. We expect that this is a general feature of the 
mechanical system defined by equations of motion of the type Eq. (2.7) with hydrodynamic 
interactions as detailed in Sec. VII. This explains the great advantage that birds in air have 
over fish in water. 

A second result of the analysis is that the mean power is equal to the mean rate of 
dissipation, as shown in Sec. VI. There is no additional energy loss related to the rate of 
change of the virtual mass. 

The mean speed and mean power can be evaluated for more sophisticated model swimmers 
or flyers by similar analysis. Elsewhere we studied longer chains with both longitudinal and 


(9.1) 
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transverse excitation in the resistive limit [2^. That analysis can be extended to the fnll 


range of viscosity, based on the equations of motion Eqs. (2.6) and (2.7). 

For assumed periodic displacements the velocity of the assembly can be derived from 
the equation of motion Eq. (3.8) by decomposition in harmonics, as demonstrated in the 
three-sphere model. The mean power is found as a collateral of the calculation. The full 
range of viscosity is covered, so that the analysis provides an interesting link between the 
flying of birds, the swimming of hsh, and the swimming of bacteria. 
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Figure captions 


Fig. 1 

Plot of the eigenvalue A+ and the absolute value |.^+| for the corresponding eigenvector 
= (1, ,^_|_) of the three-sphere model with d = 5a as functions of the dimensionless viscosity 
V*- 


Fig. 2 

Snapshot of the spheres and their velocities in the rest frame at t = 0 for the motion 
given by Eq. (8.2) with e = 3, d = 5a and rj^ = 0.01. 


Fig. 3 

Plot of the elliptical orbit in the rir 2 plane corresponding to Eq. (8.2) with e = 0.1, d = 
5a and ? 7 * = 0.01 (solid curve). We also plot the elliptical orbit for the high viscosity limit 
(dashed curve). The two curves cannot be distinguished on the scale of the hgure. 


Fig. 4 

Plot of the reduced mean swimming velocity U/{e‘^uja) for d = 5a and p* = 0.01 as a 
function of the amplitude e as calculated from Eq. (3.8) for displacements given by Eq. 
(8.3). 


Fig. 5 

Parametric plot of the reduced mean swimming power V/{e^rfuj'^a?) vs. the reduced mean 
swimming velocity U /(e^caa) for d = 5a, = 0.01 and 0 < e < 3. 

Fig. 6 

As in Fig. 4 for the efficiency Et = rjUJO^UI'D. 


Fig. 7 

Plot of the impetus X{t) (solid curve), the resistive contribution to the impetus (short 
dashes) and of the center velocity U{t) (long dashes) as functions of time during a period 
for the three-sphere swimmer for d = 5a, p* = 0.01 and amplitude factor e = 3. 
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Fig. 8 


Plot of the absolute values of the Fourier coefficients /„ of harmonics of frequency ncu, 
normalized to /o = 1, of the center velocity U{t) of the three-sphere swimmer for d = 5a 
and ? 7 * = 0.01 with amplitude factor e = 3. 


Fig. 9 

Plot of the positions of the three spheres found by numerical integration of the equations 
of motion Eq. (2.7) for d = 5a, ? 7 * = 0.01, a = 10, e = 1.5 for hfty periods of time. 


Fig. 10 

Plot of the mean value of the kinetic energy for successive periods k = 1,...,50 corre¬ 
sponding to Fig. 9. 


Fig. 11 

Plot of the orbit during the last period of Fig.9 (solid curve) compared with the elliptical 
orbit given by Eq. (8.2) (dashed curve). 
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